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Abstract 

We derive the averaged equations describing a series array of Josephson 
junctions shunted by a parallel inductor-resistor-capacitor load. We assume 
that the junctions have negligable capacitance (/? = 0), and derive averaged 
equations which turn out to be completely tractable: in particular the 
stability of both in-phase and splay states depends on a single parameter, 
5. We find an explicit expression for 6 in terms of the load parameters 
and the bias current. We recover (and refine) a common claim found in 
the technical literature, that the in-phase state is stable for inductive loads 
and unstable for capacitive loads. 
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1 Introduction 



Josephson junction arrays are perhaps the most widely studied class of coupled 
nonlinear oscillator systems. This stems in large part from their relevance in 
a number of applications, including their use as voltage standards |I| and their 
potential as sub-millimeter wave generators^] and parametric amplifiers 0. They 
also serve as prime examples of nonlinear dynamical systems with many degrees 
of freedom. Particularly good progress has been made for a subclass of this 
category, namely globally coupled oscillators. Examples of this type — where 
each oscillator is coupled with equal strength to all others — arise not only in 
the context of electrical circuits, but in the fields of laser physics and classical 
mechanics as well. 

Recent theoretical work has shown that some Josephson arrays have remark- 
able dynamical properties. The most striking discovery was made by Watanabe 
and Strogatz |Q for the class of arrays depicted in fig. |l], namely one-dimensional 
series arrays of N identical zero-capacitance junctions, driven by a constant cur- 
rent and shunted by a parallel load. Using a clever change of coordinates they 
showed that the differential equations admit — 3 independent constants of mo- 
tion, for any > 3. Furthermore, they found a rigorous reduction of the problem 
to a 5-dimensional system of differential equations, independent of N. 

The same technique allowed Watanabe and Strogatz to completely analyze 
the dynamics of the A^-oscillator system 

ipi = l + — Y,cos{ipj -ipi~6) , i = l,...,N. (1) 
i=i 

They observed in particular that the central issue — whether the attracting 
dynamics is the in-phase (i.e. synchronized) oscillation or an incoherent state — 
depends only on the sign of k sin 6. 

The purpose of the present paper is to derive Eq. (|I|) as the averaged version 
of the Josephson junction array shown in figure |l|, and to obtain an explicit 
expression for the key parameters k and 6. Our approach follows closely that 
of reference 0, which treated the special case of a pure resistive load. Starting 
from the full circuit equations, we apply a first-order averaging method which 
is valid in the weakly-coupled limit, but holds for a general inductor- resistor- 
capacitor (LRC) load. We also add some new observations about the behavior 
of the averaged system. A nice feature of the present analysis is that the results 
admit a direct physical interpretation: The combined current of the Josephson 
junction oscillators acts as a periodic driving voltage for the LRC circuit, and the 
m-phase oscillation is stable when the Josephson junction frequency is larger than 
the resonant frequency of the LRC circuit. Conversely, if the junction frequency 
is smaller than the resonant frequency of the load, then the manifold of incoherent 
states is stable. Thus we recover (and refine) an oft quoted piece of conventional 
wisdom found in the Josephson array literature, that the in-phase state is stable 
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for inductive loads, but unstable for capacitive loads. Furthermore, we observe 
that the way to extract the most energy from the Josephson junctions is to tune 
their frequency to match the resonant frequency of the LRC circuit. 

Strogatz and MiroUo [^j computed analytically the stablility of the splay-phase 
states (the most symmetric of the incoherent states) of the unaveraged system in 
the — > oo limit. Surprisingly, their results showed excellent agreement with nu- 
merical calculations even for N as small as four. Though Strogatz and MiroUo 
did not assume weak coupling, their stability results do not easily allow a physical 
interpretation. We show that our results derived from the averaged system (for 
any N) agrees with their results (for N ^ oo) in the weak coupling limit. The 
weak coupling limit also provides us with a simple physical interpretation of these 
results. 



2 Derivation of the averaged equations 

Consider the array depicted in fig. |I|. The goal of this section is to show that, in 
the limit of large shunt impedence, the circuit dynamics is governed by differential 
equations of the form ([^), and to evaluate k and S in terms of the physical 
parameters of the system. 

Our starting point is the Kirchoff equations for the circuit. We assume that the 
junctions are identical and have negligible capacitance (/5 = in the common 
notation). The governing circuit equations are 



2eRj 



+ I^sm(l)k + Q = h (2) 



LQ + RQ + ^ = -Y:<P, (3) 

where k=l, 2, . . . , N. Here, 0fc is the quantum phase difference across the 
k*^^ Josephson junction, Rj is the junction resistance, Ic is the junction critical 
current, Q is the charge on the load capacitor. If, is the applied bias current, L, 
R, and C are the load inductance, resistance, and capacitance, respectively; h is 
Planck's constant divided by 27r, e is the electron charge, and the overdot denotes 
differentiation with respect to time t. Substitution of Eq. (|) into Eq. (||) yields 

LQ + {R + NRj)Q + ^ = NRjh - Rjl, sin 0,- (4) 

j=i 

It is convenient to shift the load variable Q by a constant 

§ - NRjh - § (5) 

so that this becomes 
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■ o ^ 

LQ + {R + NRj)Q + ^ = -Rjh E sm 0i (6) 

i=i 

In order to compare arrays having different numbers of junctions, it is natural 
to define scaled load parameters: 

/ = L/N; r = R/N; c = NC (7) 
Introducing the dimensionless time r and charge q defined by 

r = jRjIct (8) 

q = lRjI^{^^\ (9) 
the circuit equations (@) and (H) become 

0fc + sin0fc + eg = a (10) 
1 ^ 

g + 7g + ujQ^q = -— ^ sin (pj (11) 

i=i 

where the overdot now denotes differentiation with respect to dimensionless time 
r, and where 

a = h/Ic (13) 
7 = (14) 

Note that 7 includes the effect of both the junction resistance and the load resis- 
tance. Thus 7 is never zero: in fact 7 > e with equality when the load resistance 
is zero. 

Up to this point, Eqs. (p!OD, (0) are merely scaled versions of the exact circuit 
equations (H), (^. To get things in a form suitable for averaging, we transform 
from the variables 0^ to natural angles ipk ||5[. The latter are "natural" in the 
sense that, in the uncoupled limit, the angular velocity 0^ is non-uniform, while 
ipk is a constant. This is accomplished by the transformation O 



"0(0) = 2 arctan 



la — 1 ((f) 7T 
^'"^2 + 4. 



(16) 
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(piip) = 2 arctan 



la + 1 



tan 



a 



TC 



Equation (^) becomes 



--i'k 



eq 



a" 



a 



sm 



(17) 



(18) 



Note that uj = i/a^ — 1 is the frequency of an uncoupled junction. It is convenient 
to rescale time to units where this frequency is unity, by taking rya^-l — > t, 
so that Eqs.(|l8D,(|ll]) become 



1 



ecuq 



a — sm 



and 



N 



uj'^q + '^ujq + ujQ'q 



T7Esin0(^i' 



N 



(19) 



(20) 



The system of equations 



, [2^) is exact, and so far we have not made an 
assumption of weak couphng. Notice that Eq.(|T^ is of the form ip^ = 1 + 0(e). 
For small e, we can obtain the drift of ipk by averaging Eq. ([T9|) over one period 
of the oscillation to obtain 



1 

1 , 

277 Jo 



2tt 



euq 



-dr. 



a 



(21) 



sin0(^fcj 

We can proceed by using the function g(r) obtained by solving Eqs.(|l9D and (^) 
with e = 0. Mathematically, the resulting (first-order) averaged system (pT]) is 
valid in the case where the Josephson junctions are weakly coupled, i.e. when 
e <^ 1. Physically, Eq.([T^) shows this occurs for sufficiently large values of the 
load inductance (per junction), and a small current flows through the load. 

To flnd the appropriate expression for g('r), we begin with the e = solution 
to Eq.(|l9D: 



'ipkir) = T + Ck 



(22) 



where the Ck are arbitrary inital conditions. From Eq. ( |T7| ) there follows the 
useful trigonometric identity 



sin (f){ip) = a — 



a' -I 



a — cos ip 

It is convenient to write this even function of ip in terms of its Fourier series 



(23) 



sin (pi^ip) = An cos^nip) 

n=0 



(24) 



where, in particular. 
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aV — 1) 



(25) 



Note that Ax is a decreasing function of a with > > —1 for a > 1. 
Combining Egs.®, (|2|)-(|2|) yields 



TV cxD 

TrEE^n cos n(r + Cj] 

j=\ n=0 



(26) 



This equation has the steady state solution 



N oo 



N 



COs(n(r + Cj) + pr, 



j=l n=0 



where 



A'^ 



(n^cu^ — uj^y + {'-fnujy 



arctan 



(27) 



(28) 



(29) 



Note that Bn and /3„ are just the well-known amplitude and phase shift response 
of a linear damped oscillator driven at frequency mo. The relative sign betweenthe 
correct branch of the inverse tangent. We will choose B^ to be positive, and since 
Ax is negative, we have < /3i < vr. 

The next step is substitute this expression for g(r) back into Eq. (pT]) . Note 
that the identity Eq. (|2^) allows us to rewrite Eq. (|2l|) as 



2tt 



a - cos(r + Ck^ 



■dr 



(30) 



Now, since g(r) is 27r-periodic, it is evident from Eq. (|30D that only the funda- 
mental Fourier component of g(r) contributes to the integral, with the result 



eBx 



N 



i^k) = 1 + sin (cj -Ck- I3x) 



(31) 



The final step is to replace the "initial values" Ck by their slowly evolving coun- 
terparts {ipki'T')), and drop the angular brackets to get the first-order averaged 
equations 



eB ^ 

^k = i + cos {ipj --ipk-s) 



(32) 



where 5 = | + /3i is given by 

sin 6 



LU — UJn 



(33) 



with f < (5 < f . This is the main result of the paper. Note that S is simply 
related to which is the phase shift of the linear LRC circuit (Eq.(p!OD) driven 
at the frequency (u) of the uncoupled Josephson Junctions (Eq. ([lT|) ). 

In terms of the parameters in the original equations, (Q) and the important 
dimensionless frequencies are: 



3 In-phase and incoherent states 

We start this section by recalling some results of references 0,0, and . Equa- 
tion (|3^) admits two types of solutons: The in-phase (or coherent) solution, where 
all of the angles are equal, and incoherent solutions, where the "center of mass" 
of the N angles ipi, when they are placed on a circle, is at the center of the circle. 
The incoherent solutons rotate rigidly and have period exactly 27r. (It is easy to 
see that the coupling terms in equation (|32|) cancel out for incoherent solutions.) 

The in-phase solution is unique: In geometric terms, this periodic orbit is a 
circle (a one- dimensional manifold) in the + 2 dimensional phase space. It is 
natural to ask "How many incoherent solutions are there?" There is a unique 
incoherent solution if iV = 2 or 3, but an infinite number of incoherent solutions 
if > 3. In fact the set of incoherent solutions is an — 2 dimensional manifold 
for any > 3, called the incoherent manifoldhj Watanabe and Strogatz ||^. The 
incoherent manifold is foliated by circles, which are the the incoherent solutions. 

We can compute the dimension of the incoherent maifold as follows: Place 
N — 2 oscillators on the circle of radius 1, so that their center of mass is not at 
the origin. (This gives the N — 2 dimensions of the incoherent manifold.) Then, 
the position of the last 2 oscillators is uniquely determined since the center of 
mass of all A^ oscillators is at the origin. Note that if the center of mass of the 
first N — 2 oscillators is farther than 2 / ( A^ — 2) from the origin, then it is not 
possible to place that last two oscillators so that the center of mass of all A^ is at 
the origin. 

A major result of reference |^ is that unaveraged systems with a "sinusoidal" 
nonlinearity, including Eqs. (^) and (^, have an incoherent manifold. In other 
words, any initial condition in the incoherent manifold is part of a periodic orbit. 
Incoherent solutions can be defined for unaveraged systems in terms of time 
delays|,|]. 

The most symmetric of the incoherent solutions, with the angles all equally 
spaced in time, is called the splay-phase solution. The stability of the in-phase 
and splay-phase solutions in equation (|3^) is easy to calculate, following reference 
(section 6). We give the stability of these periodic solutions in terms of the 
Floquet exponents, which are analogous to the eigenvalues of the linearization 
about a fixed point. A given periodic orbit has as many exponents as there are 
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phase space dimensions; if any of the exponents has positive real part, then a 
typical perturbation will grow exponentially, and the periodic orbit is unstable. 
All periodic orbits have at least one zero Floquet exponent corresponding to a 
perturbation along the orbit. 

The in-phase solution has a single Floquet exponent equal to zero, and — 1 
exponents equal to 

\n-phase = — ^ (35) 

where S is given by Eq.(^) and Bi is 



from Eqs. (|25D and (|28|) . The splay-phase solution has N — 2 Floquet exponents 
equal to zero and a complex conjugate pair 

€B 

Kplay = ^ (si^^ (5 ± i COS 5) (37) 

We see the crucial role of 6. If sin 6 > then the in- phase solution is stable and 
the splay-phase is unstable. If sin 5 < then in-phase solution is unstable and the 
spay-phase solution is neutrally stable. Watanabe and Strogatz M showed that 



every incoherent solution of the averaged equation (|32D has the same stability 
as the splay-phase solution. (The sign of each Floquet exponent is the same, 
although the magnitude varies.) An incoherent solution is neutrally stable to the 
A^ — 3 perturbations which leaves it in the incoherent manifold. If sin 6 < then 
perturbations normal to the incoherent manifold decay away, and we say that the 
incoherent manifold is stable. The incoherent manifold is stable exactly when the 
splay-phase state is neutrally stable. 

From Eq.(PBD we see that the in-phase solution is stable when u > uo, and the 
incoherent manifold is stable when u < uq. There is no bistability in the averaged 
equations. We note that numerical simulations of the unaveraged array equations 
show bistability in an LC-shunted Josephson array []T0|, |TT] (i.e. Equations (0) 



and (^ with R = 0). Therefore, the averaged equations have somewhat different 
dynamics than the original system, though both have an incoherent manifold of 
solutions. 

Strogatz and Mirollo computed the Floquet exponents of the splay-state of 
Eqs.(0) and (|^) in the limit A^ — > oo. They found that all but four of the Floquet 
exponents are zero, (the fact that Watantabe and Strogatz later proved is true 
for finite A^), while the other four are the eigenvalues of a 4 x 4 matrix. In our 
notation, the 4*^* order characteristic equation is 

+ 7A^ + {cuq^ + uj^)\^ + [ecu(Vcj2 + 1 - u) + -fuj^]\ + cuqW = (38) 

These eigenvalues derived for the case N ^ oo are in excellent agreement with 
numerical calculations Il7| even for A^ = 4, and it is conceivable that the result 
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is exactly true for any N. In any event, it makes sense to compare this formula 
with the eigenvalues derived from the averaged system, as we now do. 
If we set e = 0, the charatcteristic equation Eg . (1381) factors 

(A2 + 7A + a7o')(A' + ^') = 0. (39) 

which has a natural physical interpretation. Namely, the first factor corresponds 
to the decay of the current in the LRC branch of the circuit, while the pure 
imaginary eigenvalues ±za; from the second factor correspond to the oscillation 
frequency of a single junction in the absence of a load. If one computes from 
Eq.(^) the order e correction to the Floquet exponents ±iuj, one finds precisely 
the result for the averaged system, namely Eq. (|H7|) . Physically, then, the averaged 
result captures to lowest order the effect of interactions between the load and the 
junctions. 

Finally, we can see how accurately the averaged equations capture the transi- 
tion where the splay states go unstable. In terms of our own parameters, equation 
(14) of reference H gives the transition curve 



UJo"^ = LO^ + -uj{V^j^+T - oo) (40) 

7 

Recalling that the averaged system has the corresponding transition at ujq = uj, 
we see there is exact agreement only for e = 0. However, for e > the transition 
curve determined from Eq.(^0[) never gets very far from the diagonal when plotted 
on the u-ujq plane, as shown in fig. ^ Note from Eqs.(|l2D and (0) that e/7 = 



Rj/{r + Rj) so that this ratio never gets too big: < e/7 < 1. Moreover, for large 
uj (i.e. the limit of large bias current /{,) the transition curve always approaches 
the line u = ujq. 



4 Discussion 

Our main result is the derivation of the averaged system ( |32[ ) from the original 
Eqs.(0J^). Watanabe and Strogatz have shown that while the original equations 
are unusually tractable owing to the existence of a great many constants of mo- 
tion, the averaged equations are completely solvable . In this paper we derived 
an explicit formula for the coupling-phase 6, which is the key parameter governing 
the stability of both the in-phase and splay states. 

An advantage of the averaged equations is that they admit a fairly direct 
physical interpretation of the main features of the array dynamics. Usually, for 
a series LRC combination one identifies two resonance frequencies, which we can 
call the natural resonance frequency uq and the shifted resonance frequency ui. 
In terms of our dimensionless units, we have 

u^ = Uo- 7V2 (41) 
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where uq and 7 were defined earlier via Eqs.(|I3]) and (|Hp, respectively. In the 
presence of a periodic driving voltage at a frequency uj, the current oscillations 
(and thus the power dissipated in the load) is greatest when u = uq, while the 
capacitor's charge oscillations are greatest when u = Ui. Of course, for "high-Q" 
circuits, the current- and charge-response curves are sharply peaked and ujq is 
very nearly equal to Ui. 

Now, we noted in the last section that the averaged Josephson array equations 
display a single dynamical transition at a; = c^o- It follows that the maximum 
power delivered to the LRC load is attained at the transition point. We also recover 
an old result found in the Josephson array literature |l2|, namely that stable in- 
phase operation of a series array requires that the load "looks inductive" : an LRC 
load is said to look inductive if its impedance has positive imaginary part, which 
is equivalent to the condition oj > ujq. We note that the averaged equations admit 
an analogous stability principle for the splay state: the splay state is stable if the 
load impedance has negative imaginary part {i.e. if the load looks capacitive). 
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Figures 



Figure 1: Circuit schematic for the shunted Josephson array. 



Figure 2: Bifurcation of the splay state as given by Eq.(^OD for e/7 = 0, 0.25, 
0.5, and 0.75, and 1. The curve with e = is the resuh for the averaged system. 
The uppermost curve, with e/7 = 1, corresponds to an LC load (i? = 0). 



12 



Wiesenfeld and Swift. "Averaged Equations...", Fig.1 





Wiesenfeld and Swift. "Averaged Equations...", Fig. 2 



